# Title: Spatial Autocorrelation Analysis for Raster Data

# Objective:
# - Evaluate the spatial dependency of raster data using Moran's I calculation through multiple methods.

# Authors:
# - Oliver Gutiérrez-Hernández: http://orcid.org/0000-0003-2580-5465. E-mail: olivergh@uma.es
# - Luis V. García: http://orcid.org/0000-0002-5514-2941. E-mail: lv.garcia@csic.es

# Source:
# - Gutiérrez-Hernández, O., & García, L.V. (2025). 
#   False Discovery Rate Estimation and Control in Remote Sensing: 
#   Reliable Statistical Significance in Spatially Dependent Gridded Data. 
#   Remote Sensing Letters, 16(5), 537–554. https://doi.org/10.1080/2150704X.2025.2478664

# References for Spatial Autocorrelation and Moran's I:
# - Moran, P. A. P. 1950. "Notes on Continuous Stochastic Phenomena". Biometrika, 37(1–2): 17–23. https://doi.org/10.1093/biomet/37.1-2.17

# Practical Applications:
# - Understanding the spatial structure of data and its dependencies.
# - Understanding spatial autocorrelation is particularly relevant because FDR methods, are typically applied under the assumption of moderate positive dependency. 
# - Assessing clustering, dispersion, or randomness in raster datasets.
# - Useful for environmental studies, geostatistics, and remote sensing applications where spatial relationships are critical.

# Input:
# - Raster file: A georeferenced TIFF file containing numerical values to analyse spatial dependency.
# - Provide the full path to the raster file. Example: "C:/data/p-values.tif".
# - Ensure the raster file meets the following criteria:
#   - Contains valid numerical data.
#   - Missing or invalid values (e.g., NA) should be clearly marked.

# Outputs:
# 1. Summary Statistics:
#    - Total number of pixels.
#    - Number of valid and invalid pixels.
# 2. Moran's I Results:
#    - Global Moran's I calculated using:
#      - `terra` package
#      - `elsa` package
#      - `geocmeans` package
#      - `raster` package
# 3. Visualisation:
#    - Optional: Spatial visualisation of Moran's I results (to be implemented if needed).

# Notes:
# - This script handles missing values by replacing them with the mean of valid pixels.
# - The results from different methods may vary slightly due to differences in implementation.
# - Ensure that dependencies (`terra`, `elsa`, `geocmeans`, and `raster`) are installed and loaded before execution.


# Start of the script


# --- Prerequisites --- # !!! ATTENTION: ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ 

install_if_missing <- function(package) {
  if (!requireNamespace(package, quietly = TRUE)) {
    install.packages(package, dependencies = TRUE)
  } else {
    cat(paste("Package", package, "is already installed.\n"))
  }
}

install_if_missing("terra")
install_if_missing("elsa")
install_if_missing("geocmeans")
install_if_missing("raster")

# Load libraries
library(terra)       # For handling raster data
library(elsa)        # For Moran's I calculation
library(geocmeans)   # For calc_moran_raster function
library(raster)      # For Moran calculation with RasterLayer

# --- Set paths and load data ---
raster_path <- "C:/data/p-values.tif"  # !!! ATTENTION: REPLACE WITH YOUR FILE PATH ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲
raster_pvalues <- rast(raster_path)    # Load the raster data using terra
raster_base_name <- tools::file_path_sans_ext(basename(raster_path))  # Extract the base name for file outputs

# Display basic raster information
cat("Total number of pixels:", ncell(raster_pvalues), "\n")
p_values <- values(raster_pvalues)

if (is.null(p_values)) {
  stop("The raster contains no valid values or is improperly loaded.")
}

valid_indices <- !is.na(p_values)
num_valid_pixels <- sum(valid_indices)
num_invalid_pixels <- sum(is.na(p_values))
cat("Number of valid pixels:", num_valid_pixels, "\n")
cat("Number of invalid pixels:", num_invalid_pixels, "\n")

# Handle NA values by replacing them with the mean of valid pixels
if (num_valid_pixels > 0) {
  p_values[is.na(p_values)] <- mean(p_values, na.rm = TRUE)
  values(raster_pvalues) <- p_values
} else {
  stop("The raster contains no valid pixels to process.")
}

# --- Moran's I Calculation ---

# 1. Using terra
if (ncell(raster_pvalues) > 1) {
  global_moran_terra <- autocor(raster_pvalues, method = "moran", global = TRUE)
  p_value_terra <- 2 * (1 - pnorm(abs(global_moran_terra)))  # Approximation of p-value
} else {
  global_moran_terra <- NA
  p_value_terra <- NA
  cat("terra: Unable to calculate Moran's I due to insufficient raster size.\n")
}

# 2. Using elsa
# Calculate Moran's I using elsa
# Convert raster to RasterLayer for elsa compatibility
if (ncell(raster_pvalues) > 1) {
  raster_for_elsa <- as(raster_pvalues, "Raster")
  global_moran_elsa <- moran(raster_for_elsa, d1 = 0, d2 = 1)  # Using elsa's moran
  p_value_elsa <- 2 * (1 - pnorm(abs(global_moran_elsa)))  # Approximation of p-value
} else {
  global_moran_elsa <- NA
  p_value_elsa <- NA
  cat("elsa: Unable to calculate Moran's I due to insufficient raster size.\n")
}

# 3. Using geocmeans
if (ncell(raster_pvalues) > 1) {
  weight_matrix <- matrix(1, nrow = 3, ncol = 3)
  global_moran_geocmeans <- calc_moran_raster(raster_pvalues, weight_matrix)
  p_value_geocmeans <- 2 * (1 - pnorm(abs(global_moran_geocmeans)))  # Approximation of p-value
} else {
  global_moran_geocmeans <- NA
  p_value_geocmeans <- NA
  cat("geocmeans: Unable to calculate Moran's I due to insufficient raster size.\n")
}

# 4. Using raster package
if (ncell(raster_pvalues) > 1) {
  r_raster <- raster(raster_pvalues)
  global_moran_raster <- Moran(r_raster)
  p_value_raster <- 2 * (1 - pnorm(abs(global_moran_raster)))  # Approximation of p-value
} else {
  global_moran_raster <- NA
  p_value_raster <- NA
  cat("raster package: Unable to calculate Moran's I due to insufficient raster size.\n")
}

# Summary of Moran's I results
cat("\n--- Summary of Moran's I ---\n")
cat("terra: ", ifelse(is.na(global_moran_terra), "NA", round(global_moran_terra, 4)), " (p-value:", ifelse(is.na(p_value_terra), "NA", round(p_value_terra, 4)), ")\n")
cat("elsa (corrected): ", ifelse(is.na(global_moran_elsa), "NA", round(global_moran_elsa, 4)), " (p-value:", ifelse(is.na(p_value_elsa), "NA", round(p_value_elsa, 4)), ")\n")
cat("geocmeans: ", ifelse(is.na(global_moran_geocmeans), "NA", round(global_moran_geocmeans, 4)), " (p-value:", ifelse(is.na(p_value_geocmeans), "NA", round(p_value_geocmeans, 4)), ")\n")
cat("raster package: ", ifelse(is.na(global_moran_raster), "NA", round(global_moran_raster, 4)), " (p-value:", ifelse(is.na(p_value_raster), "NA", round(p_value_raster, 4)), ")\n")
